load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

procedure good_way()
begin

    sigm = (/0.00133607736846712, 0.00513899915929914, 0.0116468036573691, \
             0.021517142227496,   0.0346602323216384,  0.0503740850440076, \
             0.067136789924738,   0.0824488125980373,  0.0975739588097264, \
             0.115367910909135,   0.136301572054213,   0.160928952769942, \
             0.189901763775095,   0.22398679286417,    0.264086091221295, \
             0.311260796214845,   0.366759384193747,   0.432050593290324, \
             0.508862340965566,   0.599227418892952,   0.695816036901329, \
             0.787020691248154,   0.866740639339174,   0.929434205785941, \
             0.970474653327663,   0.992539942855718/)

    sigi = (/0, \
             0.00267215473693425, 0.00760584358166403, 0.0156877637330741, \
             0.0273465207219179,  0.0419739439213588,  0.0587742261666565, \
             0.0754993536828195,  0.0893982715132551,  0.105749646106198, \
             0.124986175712072,   0.147616968396354,   0.174240937143529, \
             0.205562590406661,   0.242410995321679,   0.285761187120911, \
             0.33676040530878,    0.396758363078715,   0.467342823501934, \
             0.550381858429198,   0.648072979356706,   0.743559094445952, \
             0.830482288050356,   0.902998990627992,   0.95586942094389, \
             0.985079885711436,   1/)

    pt = 2.194

    hyam = 1-sigm
    hybm = sigm

    hyai = 1-sigi
    hybi = sigi

    P0 = pt

    system("rm -rf hybpara.26.nc")
    f = addfile("hybpara.26.nc", "c")

    setfileoption(f, "DefineMode", True)

    f_atts = True
    f_atts@comments = "This hybrid sigma-pressure coordinate is equivalent to sigma coordinate"
    fileattdef(f, f_atts)

    dim_names = (/"lev","ilev"/)
    dim_sizes = (/26,27/)
    dim_unlim = (/False,False/)
    filedimdef(f, dim_names, dim_sizes, dim_unlim)

    filevardef(f, "P0", "float", "ncl_scalar")
    filevardef(f, "hyam", "float", "lev")
    filevardef(f, "hybm", "float", "lev")
    filevardef(f, "hyai", "float", "ilev")
    filevardef(f, "hybi", "float", "ilev")

    f->P0 = (/P0/)
    f->hyam = (/hyam/)
    f->hybm = (/hybm/)
    f->hyai = (/hyai/)
    f->hybi = (/hybi/)

end

procedure bad_way()
begin

    hyam = (/0.00354463800000001, 0.00738881350000001, 0.013967214, \
             0.023944625, 0.0372302900000001, 0.0531146050000002, \
             0.0700591500000003, 0.0779125700000003, 0.0766070100000003, \
             0.0750710850000003, 0.0732641500000002, 0.071138385, \
             0.0686375349999999, 0.065695415, 0.0622341550000001, \
             0.0581621650000002, 0.0533716800000001, 0.0477359250000001, \
             0.041105755, 0.0333057, 0.02496844, 0.01709591, \
             0.01021471, 0.00480317500000001, 0.00126068, 0/)  

    hybm = (/0, 0, 0, 0, 0, 0, 0, 0.00752654500000002, 0.023907685, \
             0.04317925, 0.0658512450000003, 0.0925236850000004, \
             0.1239024, 0.16081785, 0.204247, 0.2553391, 0.315446300000001, \
             0.386159300000001, 0.469349500000002, 0.567218500000003, \
             0.671827850000003, 0.770606150000003, 0.856946050000001, \
             0.924845700000002, 0.969294150000001, 0.9925561/)

    hyai = (/0.00219406700000001, 0.00489520900000001, 0.009882418, \
             0.01805201, 0.02983724, 0.0446233400000002, 0.0616058700000002, \
             0.0785124300000004, 0.0773127100000002, 0.0759013100000003, \
             0.0742408600000002, 0.0722874400000002, 0.0699893299999998, \
             0.06728574, 0.06410509, 0.0603632200000002, 0.0559611100000001, \
             0.0507822500000001, 0.0446896000000001, 0.0375219099999999, \
             0.0290894900000001, 0.02084739, 0.01334443, 0.00708499000000001, \
             0.00252136, 0, 0/)

    hybi = (/0, 0, 0, 0, 0, 0, 0, 0, 0.01505309, 0.03276228, 0.05359622, \
             0.0781062700000006, 0.1069411, 0.140863700000001, 0.180772, \
             0.227722, 0.282956200000001, 0.347936400000002, 0.4243822, \
             0.514316800000003, 0.620120200000002, 0.723535500000004, \
             0.817676800000001, 0.896215300000001, 0.953476100000003, 0.9851122, 1/)

    p0 = 1.0e5
    ps0 = 101304.196538893
    pt = 219.406700000001

    sigm = new(dimsizes(hyam), float)
    sigi = new(dimsizes(hyai), float)

    do k = 0, dimsizes(sigm)-1
        sigm(k) = (hyam(k)*p0+hybm(k)*ps0-pt)/(ps0-pt)
    end do
    do k = 0, dimsizes(sigi)-1
        sigi(k) = (hyai(k)*p0+hybi(k)*ps0-pt)/(ps0-pt)
    end do
    print(sigm)

    system("rm -rf hybpara.26.nc")
    f = addfile("hybpara.26.nc", "c")

    setfileoption(f, "DefineMode", True)

    f_atts = True
    f_atts@comments = "This hybrid sigma-pressure coordinate is NOT equivalent to sigma coordinate"
    fileattdef(f, f_atts)

    dim_names = (/"lev","ilev"/)
    dim_sizes = (/26,27/)
    dim_unlim = (/False,False/)
    filedimdef(f, dim_names, dim_sizes, dim_unlim)

    filevardef(f, "P0", "float", "ncl_scalar")
    filevardef(f, "hyam", "float", "lev")
    filevardef(f, "hybm", "float", "lev")
    filevardef(f, "hyai", "float", "ilev")
    filevardef(f, "hybi", "float", "ilev")

    f->P0 = (/p0/)
    f->hyam = (/hyam/)
    f->hybm = (/hybm/)
    f->hyai = (/hyai/)
    f->hybi = (/hybi/)

end

begin

    bad_way()

end
